# ============================================================================
# Replication Code for "Alphabet Soup: Randomized Ballot Order and
# the Representation of Marginalized Candidates"
# Sean Freeder, Justin de Benedictis-Kessner, and Rachel Bernhard
# ============================================================================

# Load required packages
library(tidyverse)
library(fixest)
library(broom)
library(scales)
library(janitor)

# Set working directory to replication folder
# setwd("path/to/replication")

# ============================================================================
# LOAD DATA ==================================================================
# ============================================================================

# Load the main dataset (assumes ceda_replication.rds is in the replication folder)
ceda_ordered <- read_rds("ceda_replication.rds")

# Define color scheme
vpurple = "#440154FF"
vyellow = "#FDE725FF"
vgreen = "#21908CFF"

# ===============================================================================
# FIGURE 1: Effect of First Ballot Position (Vote Share and Win Probability) ====
# ===============================================================================

# Run regression models with progressively more FEs - Vote Share
ceda_fit_allcands <- feols(new_voteshare ~ as.numeric(ballot_order==1) | 0,
                           cluster= ~ new_raceid,
                           data=filter(ceda_ordered, n_cands>1))
ceda_fit_allcands_ncands_fe <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands,
                                      cluster= ~new_raceid,
                                      data=filter(ceda_ordered, n_cands>1))
ceda_fit_allcands_yearcounty_fe <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co,
                                          cluster= ~new_raceid,
                                          data=filter(ceda_ordered, n_cands>1))
ceda_fit_allcands_yearcountyoffice_fe <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                                                cluster= ~new_raceid,
                                                data=filter(ceda_ordered, n_cands>1))
ceda_fit_allcands_yearcountyoffice_noninconly_fe <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                                                          cluster= ~new_raceid,
                                                          data=filter(ceda_ordered, n_cands>1 & incumb=="N"))

# Run regression models with progressively more FEs - Win Probability
ceda_fit_win_allcands <- feols(win ~ as.numeric(ballot_order==1) | 0,
                               cluster= ~ new_raceid,
                               data=filter(ceda_ordered, n_cands>1))
ceda_fit_win_allcands_ncands_fe <- feols(win ~ as.numeric(ballot_order==1) | n_cands,
                                         cluster= ~new_raceid,
                                         data=filter(ceda_ordered, n_cands>1))
ceda_fit_win_allcands_yearcounty_fe <- feols(win ~ as.numeric(ballot_order==1) | n_cands + year + co,
                                             cluster= ~new_raceid,
                                             data=filter(ceda_ordered, n_cands>1))
ceda_fit_win_allcands_yearcountyoffice_fe <- feols(win ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                                                   cluster= ~new_raceid,
                                                   data=filter(ceda_ordered, n_cands>1))
ceda_fit_win_allcands_yearcountyoffice_noninconly_fe <- feols(win ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                                                              cluster= ~new_raceid,
                                                              data=filter(ceda_ordered, n_cands>1 & incumb=="N"))

# Combine vote share results
effect_tab_vs <- bind_rows(
  tidy(ceda_fit_allcands) %>% mutate(model="Basic Model"),
  tidy(ceda_fit_allcands_ncands_fe) %>% mutate(model="    ... + # of candidates FEs"),
  tidy(ceda_fit_allcands_yearcounty_fe) %>% mutate(model="    ... + year and county FEs"),
  tidy(ceda_fit_allcands_yearcountyoffice_fe) %>% mutate(model="    ... + office FEs"),
  tidy(ceda_fit_allcands_yearcountyoffice_noninconly_fe) %>% mutate(model="    ... + non-incumbents only")
) %>%
  filter(term =="as.numeric(ballot_order == 1)") %>%
  mutate(model = factor(model, levels=c("    ... + non-incumbents only","    ... + office FEs",
                                        "    ... + year and county FEs","    ... + # of candidates FEs","Basic Model"), ordered = T),
         outcome="Voteshare")

# Combine win probability results
effect_tab_win <- bind_rows(
  tidy(ceda_fit_win_allcands) %>% mutate(model="Basic Model"),
  tidy(ceda_fit_win_allcands_ncands_fe) %>% mutate(model="    ... + # of candidates FEs"),
  tidy(ceda_fit_win_allcands_yearcounty_fe) %>% mutate(model="    ... + year and county FEs"),
  tidy(ceda_fit_win_allcands_yearcountyoffice_fe) %>% mutate(model="    ... + office FEs"),
  tidy(ceda_fit_win_allcands_yearcountyoffice_noninconly_fe) %>% mutate(model="    ... + non-incumbents only")
) %>%
  filter(term =="as.numeric(ballot_order == 1)") %>%
  mutate(model = factor(model, levels=c("    ... + non-incumbents only","    ... + office FEs",
                                        "    ... + year and county FEs","    ... + # of candidates FEs","Basic Model"), ordered = T),
         outcome="Probability of Winning")

effect_tab_vs_and_win <- bind_rows(effect_tab_vs, effect_tab_win) %>%
  mutate(outcome = factor(outcome, levels=c("Probability of Winning","Voteshare"), ordered=T))

# Create Figure 1
(fig1 <- ggplot(effect_tab_vs_and_win) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=model, y=estimate, group=model), size=2) +
  geom_errorbar(aes(y=estimate, x=model, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=model), width=0) +
  geom_errorbar(aes(y=estimate, x=model, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=model), lwd=1, width=0) +
  geom_text(aes(x=model, y=estimate+0.005, group=model, label=paste0(round(estimate,3)*100,"%")), nudge_x=0.35, family="Arial Narrow") +
  facet_wrap(~outcome, nrow=2, as.table = F) +
  scale_y_continuous("Effect of being first on outcome", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"),
        axis.text.y = element_text(hjust=0),
        strip.background = element_rect(fill="grey"))
)

ggsave(fig1, filename = "figures/figure1_effect_first_fes.png", height=4, width=5)

# ============================================================================
# FIGURE 2: Effect of Ballot Order by Race and Gender (Non-Incumbents) =======
# ============================================================================

# Run subgroup regressions
ceda_fit_noninconly <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                              data=filter(ceda_ordered, n_cands>1 & incumb=="N"))
ceda_fit_noninconly_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                    data=filter(ceda_ordered, n_cands>1 & incumb=="N" & female==1))
ceda_fit_noninconly_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                  data=filter(ceda_ordered, n_cands>1 & incumb=="N" & female==0))
ceda_fit_noninconly_whites <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                     data=filter(ceda_ordered, n_cands>1 & incumb=="N" & white==1))
ceda_fit_noninconly_nonwhites <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                        data=filter(ceda_ordered, n_cands>1 & incumb=="N" & white==0))
ceda_fit_noninconly_white_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                          data=filter(ceda_ordered, n_cands>1 & incumb=="N" & white==1 & female==1))
ceda_fit_noninconly_nonwhite_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                             data=filter(ceda_ordered, n_cands>1 & incumb=="N" & white==0 & female==1))
ceda_fit_noninconly_white_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                        data=filter(ceda_ordered, n_cands>1 & incumb=="N" & white==1 & female==0))
ceda_fit_noninconly_nonwhite_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                           data=filter(ceda_ordered, n_cands>1 & incumb=="N" & white==0 & female==0))

effect_tab_fig2 <- bind_rows(
  tidy(ceda_fit_noninconly) %>% mutate(race="All", gender="All"),
  tidy(ceda_fit_noninconly_women) %>% mutate(race="All", gender="Women"),
  tidy(ceda_fit_noninconly_men) %>% mutate(race="All", gender="Men"),
  tidy(ceda_fit_noninconly_whites) %>% mutate(race="White", gender="All"),
  tidy(ceda_fit_noninconly_nonwhites) %>% mutate(race="Non-white", gender="All"),
  tidy(ceda_fit_noninconly_white_women) %>% mutate(race="White", gender="Women"),
  tidy(ceda_fit_noninconly_nonwhite_women) %>% mutate(race="Non-white", gender="Women"),
  tidy(ceda_fit_noninconly_white_men) %>% mutate(race="White", gender="Men"),
  tidy(ceda_fit_noninconly_nonwhite_men) %>% mutate(race="Non-white", gender="Men")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T))

(fig2 <- ggplot(effect_tab_fig2) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=gender, y=estimate, group=race, shape=race), size=2, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=race),
                width=0, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=race),
                lwd=1, width=0, position=position_dodge(width=0.4)) +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  scale_shape_manual("Candidate Race", breaks=c("All","White","Non-white"), values=c(3,1,16)) +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0))
)

ggsave(fig2, filename = "figures/figure2_effect_racegender_nonincumbents.png", height=3, width=5)

# ============================================================================
# FIGURE 3: Effect of Ballot Order by Election Timing ========================
# ============================================================================

ceda_fit_cands_pres <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                              cluster= ~new_raceid,
                              data=filter(ceda_ordered, n_cands>1 & incumb=="N" & electiming3=="Presidential"))
ceda_fit_cands_mid <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                             cluster= ~new_raceid,
                             data=filter(ceda_ordered, n_cands>1 & incumb=="N" & electiming3=="Midterm"))
ceda_fit_cands_off <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2,
                             cluster= ~new_raceid,
                             data=filter(ceda_ordered, n_cands>1 & incumb=="N" & electiming3=="Off-cycle"))

coefs_for_plot_bytiming <- bind_rows(
  tidy(ceda_fit_cands_pres),
  tidy(ceda_fit_cands_mid),
  tidy(ceda_fit_cands_off)
) %>%
  filter(term=="as.numeric(ballot_order == 1)") %>%
  mutate(timing = rep(c("Presidential","Midterm","Off-cycle"), each=1),
         timing = factor(timing, levels=c("Presidential","Midterm","Off-cycle"), ordered=T))

(fig3 <- ggplot(coefs_for_plot_bytiming) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=timing, y=estimate, group=term), position=position_dodge(width=0.5), size=2) +
  geom_errorbar(aes(y=estimate, x=timing, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=term),
                width=0, position=position_dodge(width=0.5)) +
  geom_errorbar(aes(y=estimate, x=timing, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=term),
                lwd=1, width=0, position=position_dodge(width=0.5)) +
  geom_text(aes(x=as.numeric(timing)+0.005, y=estimate+0.001, group=term, label=paste0(round(estimate,3)*100,"%")),
            nudge_x=0.2, family="Arial Narrow") +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"),
        axis.title.y = element_text(margin = margin(r = 20)),
        plot.margin = margin(.5, .5, 0, .5, "cm"))
)

ggsave(fig3, filename = "figures/figure3_effect_bytiming.png", height=2.5, width=6)

# ============================================================================
# FIGURE 4: Effect of Ballot Order by Race, Gender, and Election Timing ======
# ============================================================================

# Presidential elections
ceda_fit_pres <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                        data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N"))
ceda_fit_pres_white <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                              data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & white==1))
ceda_fit_pres_nonwhite <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                 data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & white==0))
ceda_fit_pres_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                              data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & female==1))
ceda_fit_pres_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & female==0))
ceda_fit_pres_white_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                    data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & white==1 & female==1))
ceda_fit_pres_nonwhite_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                       data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & white==0 & female==1))
ceda_fit_pres_white_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                  data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & white==1 & female==0))
ceda_fit_pres_nonwhite_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                     data=filter(ceda_ordered, electiming3=="Presidential" & n_cands>1 & incumb=="N" & white==0 & female==0))

effect_tab_pres <- bind_rows(
  tidy(ceda_fit_pres) %>% mutate(race="All", gender="All"),
  tidy(ceda_fit_pres_white) %>% mutate(race="White", gender="All"),
  tidy(ceda_fit_pres_nonwhite) %>% mutate(race="Non-white", gender="All"),
  tidy(ceda_fit_pres_women) %>% mutate(race="All", gender="Women"),
  tidy(ceda_fit_pres_men) %>% mutate(race="All", gender="Men"),
  tidy(ceda_fit_pres_white_women) %>% mutate(race="White", gender="Women"),
  tidy(ceda_fit_pres_nonwhite_women) %>% mutate(race="Non-white", gender="Women"),
  tidy(ceda_fit_pres_white_men) %>% mutate(race="White", gender="Men"),
  tidy(ceda_fit_pres_nonwhite_men) %>% mutate(race="Non-white", gender="Men")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T),
         plotorder = as.numeric(gender),
         timing = "Presidential")

# Midterm elections
ceda_fit_midterm <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                           data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N"))
ceda_fit_midterm_white <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                 data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & white==1))
ceda_fit_midterm_nonwhite <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                    data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & white==0))
ceda_fit_midterm_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                 data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & female==1))
ceda_fit_midterm_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                               data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & female==0))
ceda_fit_midterm_white_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                       data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & white==1 & female==1))
ceda_fit_midterm_nonwhite_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                          data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & white==0 & female==1))
ceda_fit_midterm_white_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                     data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & white==1 & female==0))
ceda_fit_midterm_nonwhite_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                        data=filter(ceda_ordered, electiming3=="Midterm" & n_cands>1 & incumb=="N" & white==0 & female==0))

effect_tab_midterm <- bind_rows(
  tidy(ceda_fit_midterm) %>% mutate(race="All", gender="All"),
  tidy(ceda_fit_midterm_white) %>% mutate(race="White", gender="All"),
  tidy(ceda_fit_midterm_nonwhite) %>% mutate(race="Non-white", gender="All"),
  tidy(ceda_fit_midterm_women) %>% mutate(race="All", gender="Women"),
  tidy(ceda_fit_midterm_men) %>% mutate(race="All", gender="Men"),
  tidy(ceda_fit_midterm_white_women) %>% mutate(race="White", gender="Women"),
  tidy(ceda_fit_midterm_nonwhite_women) %>% mutate(race="Non-white", gender="Women"),
  tidy(ceda_fit_midterm_white_men) %>% mutate(race="White", gender="Men"),
  tidy(ceda_fit_midterm_nonwhite_men) %>% mutate(race="Non-white", gender="Men")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T),
         plotorder = as.numeric(gender),
         timing = "Midterm")

# Off-cycle elections
ceda_fit_offcycle <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N"))
ceda_fit_offcycle_white <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                  data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & white==1))
ceda_fit_offcycle_nonwhite <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                     data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & white==0))
ceda_fit_offcycle_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                  data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & female==1))
ceda_fit_offcycle_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & female==0))
ceda_fit_offcycle_white_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                        data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & white==1 & female==1))
ceda_fit_offcycle_nonwhite_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                           data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & white==0 & female==1))
ceda_fit_offcycle_white_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                      data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & white==1 & female==0))
ceda_fit_offcycle_nonwhite_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                         data=filter(ceda_ordered, electiming3=="Off-cycle" & n_cands>1 & incumb=="N" & white==0 & female==0))

effect_tab_offcycle <- bind_rows(
  tidy(ceda_fit_offcycle) %>% mutate(race="All", gender="All"),
  tidy(ceda_fit_offcycle_white) %>% mutate(race="White", gender="All"),
  tidy(ceda_fit_offcycle_nonwhite) %>% mutate(race="Non-white", gender="All"),
  tidy(ceda_fit_offcycle_women) %>% mutate(race="All", gender="Women"),
  tidy(ceda_fit_offcycle_men) %>% mutate(race="All", gender="Men"),
  tidy(ceda_fit_offcycle_white_women) %>% mutate(race="White", gender="Women"),
  tidy(ceda_fit_offcycle_nonwhite_women) %>% mutate(race="Non-white", gender="Women"),
  tidy(ceda_fit_offcycle_white_men) %>% mutate(race="White", gender="Men"),
  tidy(ceda_fit_offcycle_nonwhite_men) %>% mutate(race="Non-white", gender="Men")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T),
         plotorder = as.numeric(gender),
         timing = "Off-cycle")

effect_tab_fig4 <- bind_rows(effect_tab_pres, effect_tab_midterm, effect_tab_offcycle) %>%
  mutate(timing = factor(timing, levels=c("Presidential","Midterm","Off-cycle"), ordered = T))

(fig4 <- ggplot(effect_tab_fig4) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=plotorder, y=estimate, group=race, shape=race), size=2, position=position_dodge(width=0.6)) +
  geom_errorbar(aes(x=plotorder, y=estimate, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=race),
                width=0, position=position_dodge(width=0.6)) +
  geom_errorbar(aes(x=plotorder, y=estimate, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=race),
                lwd=1, width=0, position=position_dodge(width=0.6)) +
  geom_text(data=filter(effect_tab_fig4, gender=="All" & timing=="Presidential"),
            aes(x=plotorder, y=estimate+qnorm(0.975)*std.error+0.001, group=race, label=paste0(race," candidates")),
            family="Arial Narrow", position=position_dodge(width=0.6), hjust=0, size=2.5) +
  facet_wrap(~timing) +
  scale_y_continuous("Effect of being first on voteshare", labels=scales::percent_format(accuracy = 1)) +
  scale_x_continuous("", breaks=c(1,2,3), limits=c(0.7,3.3), labels=c("Men","Women","All"), position = "bottom") +
  scale_shape_manual("Candidate Race", breaks=c("All","White","Non-white"), values=c(3,1,16)) +
  coord_flip(ylim=c(-0.01,0.052)) +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0),
        strip.background = element_rect(fill="grey"), legend.position="none")
)

ggsave(fig4, filename = "figures/figure4_effect_racegender_bytiming.png", height=3, width=8)

# ============================================================================
# FIGURE 5: Effect of Ballot Order by Race and Gender (Incumbents) ===========
# ============================================================================

ceda_fit_inconly <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                           data=filter(ceda_ordered, n_cands>1 & incumb=="Y"))
ceda_fit_inconly_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                 data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & female==1))
ceda_fit_inconly_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                               data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & female==0))
ceda_fit_inconly_whites <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                  data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & white==1))
ceda_fit_inconly_nonwhites <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                     data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & white==0))
ceda_fit_inconly_white_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                       data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & white==1 & female==1))
ceda_fit_inconly_nonwhite_women <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                          data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & white==0 & female==1))
ceda_fit_inconly_white_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                     data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & white==1 & female==0))
ceda_fit_inconly_nonwhite_men <- feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                                        data=filter(ceda_ordered, n_cands>1 & incumb=="Y" & white==0 & female==0))

effect_tab_fig5 <- bind_rows(
  tidy(ceda_fit_inconly) %>% mutate(race="All", gender="All"),
  tidy(ceda_fit_inconly_women) %>% mutate(race="All", gender="Women"),
  tidy(ceda_fit_inconly_men) %>% mutate(race="All", gender="Men"),
  tidy(ceda_fit_inconly_whites) %>% mutate(race="White", gender="All"),
  tidy(ceda_fit_inconly_nonwhites) %>% mutate(race="Non-white", gender="All"),
  tidy(ceda_fit_inconly_white_women) %>% mutate(race="White", gender="Women"),
  tidy(ceda_fit_inconly_nonwhite_women) %>% mutate(race="Non-white", gender="Women"),
  tidy(ceda_fit_inconly_white_men) %>% mutate(race="White", gender="Men"),
  tidy(ceda_fit_inconly_nonwhite_men) %>% mutate(race="Non-white", gender="Men")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T))

(fig5 <- ggplot(effect_tab_fig5) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=gender, y=estimate, group=race, shape=race), size=2, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=race),
                width=0, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=race),
                lwd=1, width=0, position=position_dodge(width=0.4)) +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  scale_shape_manual("Candidate Race", breaks=c("All","White","Non-white"), values=c(3,1,16)) +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0))
)

ggsave(fig5, filename = "figures/figure5_effect_racegender_incumbents.png", height=3, width=5)

# ============================================================================
# FIGURE 6: Ballot Order, Candidate Race, and Partisanship (Bar Plots) =======
# ============================================================================

tab_for_hist <- ceda_ordered %>%
  mutate(above_avg_demshare_pres = ifelse(demshare_pres > median(ceda_ordered$demshare_pres, na.rm=T),
                                          "Above Median Dem. Presidential Voteshare",
                                          "Below Median Dem. Presidential Voteshare")) %>%
  filter(n_cands==2) %>%
  group_by(above_avg_demshare_pres, white, first) %>%
  summarize(avg_voteshare = mean(new_voteshare, na.rm=T)) %>%
  na.omit() %>%
  mutate(first_cat = as.factor(case_when(first==1 ~ "First", first==0 ~ "Not first")),
         race_cat = as.factor(case_when(white==1 ~ "White", white==0 ~ "Non-white")),
         above_avg_demshare_pres = factor(above_avg_demshare_pres,
                                          levels=c("Below Median Dem. Presidential Voteshare",
                                                   "Above Median Dem. Presidential Voteshare"), ordered=T))

(fig6 <- ggplot(tab_for_hist) +
  geom_histogram(aes(x=first_cat, y=avg_voteshare, fill=race_cat, group=race_cat), stat="identity",
                 size=2, position=position_dodge(width=1)) +
  geom_text(aes(x=first_cat, y=avg_voteshare+0.005, group=race_cat, label=paste0(round(avg_voteshare,2)*100,"%")),
            family="Arial Narrow", position=position_dodge(width=1)) +
  facet_wrap(~above_avg_demshare_pres, as.table = T) +
  scale_y_continuous("Average voteshare", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("Ballot position") +
  scale_fill_manual("Candidate race", breaks=c("Non-white","White"), values=c(vpurple, vgreen)) +
  coord_cartesian(ylim=c(0.4,0.55)) +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0), legend.position = "bottom",
        strip.background = element_rect(fill="grey"), strip.text = element_text(size=13))
)

ggsave(fig6, filename = "figures/figure6_barplot_bypvote.png", height=5, width=7)

# ============================================================================
# FIGURE 7: Effect by Race/Gender and Democratic Presidential Vote Share =====
# ============================================================================

# Helper function to run models by demshare category
run_demshare_models <- function(demshare_cat_val) {
  list(
    all = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val)),
    women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                  data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & female==1)),
    men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & female==0)),
    whites = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                   data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & white==1)),
    nonwhites = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                      data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & white==0)),
    white_women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                        data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & white==1 & female==1)),
    nonwhite_women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                           data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & white==0 & female==1)),
    white_men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                      data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & white==1 & female==0)),
    nonwhite_men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & demshare_pres_cat == demshare_cat_val & white==0 & female==0))
  )
}

combine_demshare_results <- function(models, pvote_label) {
  bind_rows(
    tidy(models$all) %>% mutate(race="All", gender="All", pvote=pvote_label),
    tidy(models$women) %>% mutate(race="All", gender="Women", pvote=pvote_label),
    tidy(models$men) %>% mutate(race="All", gender="Men", pvote=pvote_label),
    tidy(models$whites) %>% mutate(race="White", gender="All", pvote=pvote_label),
    tidy(models$nonwhites) %>% mutate(race="Non-white", gender="All", pvote=pvote_label),
    tidy(models$white_women) %>% mutate(race="White", gender="Women", pvote=pvote_label),
    tidy(models$nonwhite_women) %>% mutate(race="Non-white", gender="Women", pvote=pvote_label),
    tidy(models$white_men) %>% mutate(race="White", gender="Men", pvote=pvote_label),
    tidy(models$nonwhite_men) %>% mutate(race="Non-white", gender="Men", pvote=pvote_label)
  )
}

models_demshare_1 <- run_demshare_models("<50%")
models_demshare_2 <- run_demshare_models("50-54%")
models_demshare_3 <- run_demshare_models("55-71%")
models_demshare_4 <- run_demshare_models(">72%")

effect_tab_fig7 <- bind_rows(
  combine_demshare_results(models_demshare_1, "<50%"),
  combine_demshare_results(models_demshare_2, "50-54%"),
  combine_demshare_results(models_demshare_3, "55-71%"),
  combine_demshare_results(models_demshare_4, ">72%")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T),
         demshare_pres_cat = factor(pvote, levels=c("<50%","50-54%","55-71%",">72%"), ordered=T))

(fig7 <- ggplot(effect_tab_fig7) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=gender, y=estimate, group=race, shape=race), size=2, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=race),
                width=0, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=race),
                lwd=1, width=0, position=position_dodge(width=0.4)) +
  facet_wrap(~demshare_pres_cat, labeller = as_labeller(function(x) paste0("Dem. Presidential Vote: ",x))) +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  scale_shape_manual("Candidate Race", breaks=c("All","White","Non-white"), values=c(3,1,16)) +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0), legend.position = c(0.92,0.2),
        strip.background = element_rect(fill="grey"))
)

ggsave(fig7, filename = "figures/figure7_effect_racegender_bypvote.png", height=5, width=7)

# ============================================================================
# FIGURE 8: Effect by Race/Gender and County % Hispanic/Latino ===============
# ============================================================================

# Helper function for county pct hisp models
run_pcthisp_models <- function(pcthisp_cat_val) {
  list(
    all = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val)),
    women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                  data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & female==1)),
    men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & female==0)),
    whites = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                   data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & white==1)),
    nonwhites = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                      data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & white==0)),
    white_women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                        data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & white==1 & female==1)),
    nonwhite_women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                           data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & white==0 & female==1)),
    white_men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                      data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & white==1 & female==0)),
    nonwhite_men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_hisp_cat == pcthisp_cat_val & white==0 & female==0))
  )
}

combine_pcthisp_results <- function(models, pcthisp_label) {
  bind_rows(
    tidy(models$all) %>% mutate(race="All", gender="All", countypct=pcthisp_label),
    tidy(models$women) %>% mutate(race="All", gender="Women", countypct=pcthisp_label),
    tidy(models$men) %>% mutate(race="All", gender="Men", countypct=pcthisp_label),
    tidy(models$whites) %>% mutate(race="White", gender="All", countypct=pcthisp_label),
    tidy(models$nonwhites) %>% mutate(race="Non-white", gender="All", countypct=pcthisp_label),
    tidy(models$white_women) %>% mutate(race="White", gender="Women", countypct=pcthisp_label),
    tidy(models$nonwhite_women) %>% mutate(race="Non-white", gender="Women", countypct=pcthisp_label),
    tidy(models$white_men) %>% mutate(race="White", gender="Men", countypct=pcthisp_label),
    tidy(models$nonwhite_men) %>% mutate(race="Non-white", gender="Men", countypct=pcthisp_label)
  )
}

models_pcthisp_1 <- run_pcthisp_models("<17%")
models_pcthisp_2 <- run_pcthisp_models("17-31%")
models_pcthisp_3 <- run_pcthisp_models("32-40%")
models_pcthisp_4 <- run_pcthisp_models(">40%")

effect_tab_fig8 <- bind_rows(
  combine_pcthisp_results(models_pcthisp_1, "<17%"),
  combine_pcthisp_results(models_pcthisp_2, "17-31%"),
  combine_pcthisp_results(models_pcthisp_3, "32-40%"),
  combine_pcthisp_results(models_pcthisp_4, ">40%")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T),
         county_pct_hisp_cat = factor(countypct, levels=c("<17%","17-31%","32-40%",">40%"), ordered=T))

(fig8 <- ggplot(effect_tab_fig8) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=gender, y=estimate, group=race, shape=race), size=2, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=race),
                width=0, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=race),
                lwd=1, width=0, position=position_dodge(width=0.4)) +
  facet_wrap(~county_pct_hisp_cat, labeller = as_labeller(function(x) paste0("County % Hispanic: ",x))) +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  scale_shape_manual("Candidate Race", breaks=c("All","White","Non-white"), values=c(3,1,16)) +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0), legend.position = c(0.92,0.2),
        strip.background = element_rect(fill="grey"))
)

ggsave(fig8, filename = "figures/figure8_effect_racegender_bypcthisp.png", height=5, width=7)

# ================================================================================
# APPENDIX FIGURE A1: Size of Ballot Position Effects by Number of Candidates ====
# ================================================================================

ceda_fit_ncands_2 <- feols(new_voteshare ~ as.numeric(ballot_order==1) | year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, n_cands==2 & incumb=="N"))
ceda_fit_ncands_3 <- feols(new_voteshare ~ as.numeric(ballot_order==1) | year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, n_cands==3 & incumb=="N"))
ceda_fit_ncands_4 <- feols(new_voteshare ~ as.numeric(ballot_order==1) | year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, n_cands==4 & incumb=="N"))
ceda_fit_ncands_5 <- feols(new_voteshare ~ as.numeric(ballot_order==1) | year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, n_cands==5 & incumb=="N"))
ceda_fit_ncands_6plus <- feols(new_voteshare ~ as.numeric(ballot_order==1) | year + co + office2, cluster= ~new_raceid,
                                data=filter(ceda_ordered, n_cands>=6 & incumb=="N"))

coefs_for_plot_byncands <- bind_rows(
  tidy(ceda_fit_ncands_2),
  tidy(ceda_fit_ncands_3),
  tidy(ceda_fit_ncands_4),
  tidy(ceda_fit_ncands_5),
  tidy(ceda_fit_ncands_6plus)
) %>%
  filter(term=="as.numeric(ballot_order == 1)") %>%
  mutate(n_cands = c("2","3","4","5","6+"),
         n_cands = factor(n_cands, levels=c("6+","5","4","3","2"), ordered=T))

(figA1 <- ggplot(coefs_for_plot_byncands) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=n_cands, y=estimate, group=term), position=position_dodge(width=0.5), size=2) +
  geom_errorbar(aes(y=estimate, x=n_cands, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=term),
                width=0, position=position_dodge(width=0.5)) +
  geom_errorbar(aes(y=estimate, x=n_cands, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=term),
                lwd=1, width=0, position=position_dodge(width=0.5)) +
  geom_text(aes(x=as.numeric(n_cands)+0.005, y=estimate+0.001, group=term, label=paste0(round(estimate,3)*100,"%")),
            nudge_x=0.2, family="Arial Narrow") +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("Number of candidates in race") +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"),
        axis.title.y = element_text(margin = margin(r = 20)),
        plot.margin = margin(.5, .5, 0, .5, "cm"))
)

ggsave(figA1, filename = "figures/figureA1_effect_byncands.png", height=3, width=6)

# ============================================================================
# APPENDIX FIGURE A2: Effect of Ballot Order by Office =======================
# ============================================================================

ceda_fit_citycouncil <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co, cluster= ~new_raceid,
                               data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="City Council"))
ceda_fit_countyleg <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                             data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="County Legislator"))
ceda_fit_school <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                          data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="School Board"))
ceda_fit_mayor <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="Mayor"))
ceda_fit_treasurer <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                             data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="Treasurer"))
ceda_fit_clerk <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="Clerk"))
ceda_fit_judge <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="Judge"))
ceda_fit_prosecutor <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                              data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="Prosecutor"))
ceda_fit_sheriff <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                           data=filter(ceda_ordered, n_cands>1 & incumb=="N" & office2=="Sheriff"))

coefs_for_plot_byoffice <- bind_rows(
  tidy(ceda_fit_citycouncil),
  tidy(ceda_fit_countyleg),
  tidy(ceda_fit_school),
  tidy(ceda_fit_treasurer),
  tidy(ceda_fit_clerk),
  tidy(ceda_fit_judge),
  tidy(ceda_fit_prosecutor),
  tidy(ceda_fit_sheriff),
  tidy(ceda_fit_mayor)
) %>%
  filter(term=="as.numeric(ballot_order == 1)") %>%
  mutate(office = c("City Council","County Legislator","School Board","Treasurer","Clerk","Judge","Prosecutor","Sheriff","Mayor"),
         office = factor(office, levels=c("School Board","City Council","County Legislator","Treasurer","Clerk","Judge","Prosecutor","Sheriff","Mayor"), ordered=T),
         office = fct_reorder(office, estimate))

(figA2 <- ggplot(coefs_for_plot_byoffice) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=office, y=estimate, group=term), position=position_dodge(width=0.5), size=2) +
  geom_errorbar(aes(y=estimate, x=office, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=term),
                width=0, position=position_dodge(width=0.5)) +
  geom_errorbar(aes(y=estimate, x=office, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=term),
                lwd=1, width=0, position=position_dodge(width=0.5)) +
  geom_text(aes(x=as.numeric(office)+0.005, y=estimate+0.001, group=term, label=paste0(round(estimate,3)*100,"%")),
            nudge_x=0.2, family="Arial Narrow") +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"),
        axis.title.y = element_text(margin = margin(r = 20)),
        plot.margin = margin(.5, .5, 0, .5, "cm"))
)

ggsave(figA2, filename = "figures/figureA2_effect_byoffice.png", height=4.5, width=6)

# ============================================================================
# APPENDIX FIGURE A4: Effect by Time Period ==================================
# ============================================================================

ceda_fit_early <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & timeperiod=="1995-2004"))
ceda_fit_mid_time <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                            data=filter(ceda_ordered, n_cands>1 & incumb=="N" & timeperiod=="2005-2014"))
ceda_fit_late <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                        data=filter(ceda_ordered, n_cands>1 & incumb=="N" & timeperiod=="2015-2021"))

coefs_for_plot_bytimeperiod <- bind_rows(
  tidy(ceda_fit_early),
  tidy(ceda_fit_mid_time),
  tidy(ceda_fit_late)
) %>%
  filter(term=="as.numeric(ballot_order == 1)") %>%
  mutate(timeperiod = c("1995-2004","2005-2014","2015-2021"),
         timeperiod = factor(timeperiod, levels=c("2015-2021","2005-2014","1995-2004"), ordered=T))

(figA4 <- ggplot(coefs_for_plot_bytimeperiod) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=timeperiod, y=estimate, group=term), position=position_dodge(width=0.5), size=2) +
  geom_errorbar(aes(y=estimate, x=timeperiod, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=term),
                width=0, position=position_dodge(width=0.5)) +
  geom_errorbar(aes(y=estimate, x=timeperiod, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=term),
                lwd=1, width=0, position=position_dodge(width=0.5)) +
  geom_text(aes(x=as.numeric(timeperiod)+0.005, y=estimate+0.001, group=term, label=paste0(round(estimate,3)*100,"%")),
            nudge_x=0.2, family="Arial Narrow") +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"),
        axis.title.y = element_text(margin = margin(r = 20)),
        plot.margin = margin(.5, .5, 0, .5, "cm"))
)

ggsave(figA4, filename = "figures/figureA4_effect_bytimeperiod.png", height=4.5, width=6)

# ============================================================================
# APPENDIX FIGURE A7: Effect by County % White ===============================
# ============================================================================

# Helper function for county pct white models
run_pctwhite_models <- function(pctwhite_cat_val) {
  list(
    all = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val)),
    women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                  data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & female==1)),
    men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & female==0)),
    whites = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                   data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & white==1)),
    nonwhites = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                      data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & white==0)),
    white_women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                        data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & white==1 & female==1)),
    nonwhite_women = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                           data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & white==0 & female==1)),
    white_men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                      data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & white==1 & female==0)),
    nonwhite_men = feols(new_voteshare ~ first | n_cands + year + co + office2, cluster= ~new_raceid,
                         data=filter(ceda_ordered, n_cands>1 & incumb=="N" & county_pct_white_cat == pctwhite_cat_val & white==0 & female==0))
  )
}

combine_pctwhite_results <- function(models, pctwhite_label) {
  bind_rows(
    tidy(models$all) %>% mutate(race="All", gender="All", countypct=pctwhite_label),
    tidy(models$women) %>% mutate(race="All", gender="Women", countypct=pctwhite_label),
    tidy(models$men) %>% mutate(race="All", gender="Men", countypct=pctwhite_label),
    tidy(models$whites) %>% mutate(race="White", gender="All", countypct=pctwhite_label),
    tidy(models$nonwhites) %>% mutate(race="Non-white", gender="All", countypct=pctwhite_label),
    tidy(models$white_women) %>% mutate(race="White", gender="Women", countypct=pctwhite_label),
    tidy(models$nonwhite_women) %>% mutate(race="Non-white", gender="Women", countypct=pctwhite_label),
    tidy(models$white_men) %>% mutate(race="White", gender="Men", countypct=pctwhite_label),
    tidy(models$nonwhite_men) %>% mutate(race="Non-white", gender="Men", countypct=pctwhite_label)
  )
}

models_pctwhite_1 <- run_pctwhite_models("<10%")
models_pctwhite_2 <- run_pctwhite_models("10-44%")
models_pctwhite_3 <- run_pctwhite_models("45-66%")
models_pctwhite_4 <- run_pctwhite_models(">66%")

effect_tab_figA7 <- bind_rows(
  combine_pctwhite_results(models_pctwhite_1, "<10%"),
  combine_pctwhite_results(models_pctwhite_2, "10-44%"),
  combine_pctwhite_results(models_pctwhite_3, "45-66%"),
  combine_pctwhite_results(models_pctwhite_4, ">66%")
) %>%
  filter(term == "first") %>%
  mutate(race = factor(race, levels=c("Non-white","White","All"), ordered = T),
         gender = factor(gender, levels=c("Men","Women","All"), ordered = T),
         county_pct_white_cat = factor(countypct, levels=c("<10%","10-44%","45-66%",">66%"), ordered=T))

(figA7 <- ggplot(effect_tab_figA7) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=gender, y=estimate, group=race, shape=race), size=2, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error, group=race),
                width=0, position=position_dodge(width=0.4)) +
  geom_errorbar(aes(x=gender, y=estimate, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error, group=race),
                lwd=1, width=0, position=position_dodge(width=0.4)) +
  facet_wrap(~county_pct_white_cat, labeller = as_labeller(function(x) paste0("County % White: ",x))) +
  scale_y_continuous("Effect of being first on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("") +
  scale_shape_manual("Candidate Race", breaks=c("All","White","Non-white"), values=c(3,1,16)) +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"), axis.text.y = element_text(hjust=0), legend.position = c(0.92,0.2),
        strip.background = element_rect(fill="grey"))
)

ggsave(figA7, filename = "figures/figureA7_effect_racegender_bycountypctwhite.png", height=5, width=7)

# ============================================================================
# APPENDIX FIGURE A8: Effect for All Ballot Positions ========================
# ============================================================================

ceda_fit_first_a <- feols(new_voteshare ~ as.numeric(ballot_order==1) | n_cands + year + co + office2, cluster= ~new_raceid,
                          data=filter(ceda_ordered, n_cands>1 & incumb=="N"))
ceda_fit_firstsecond_b <- feols(new_voteshare ~ as.numeric(ballot_order==1) + as.numeric(ballot_order==2) | n_cands + year + co + office2,
                                cluster= ~new_raceid, data=filter(ceda_ordered, n_cands>2 & incumb=="N"))
ceda_fit_firstsecondthird_c <- feols(new_voteshare ~ as.numeric(ballot_order==1) + as.numeric(ballot_order==2) + as.numeric(ballot_order==3) | n_cands + year + co + office2,
                                     cluster= ~new_raceid, data=filter(ceda_ordered, n_cands>3 & incumb=="N"))
ceda_fit_firstsecondthirdfourth_d <- feols(new_voteshare ~ as.numeric(ballot_order==1) + as.numeric(ballot_order==2) + as.numeric(ballot_order==3) + as.numeric(ballot_order==4) | n_cands + year + co + office2,
                                           cluster= ~new_raceid, data=filter(ceda_ordered, n_cands>4 & incumb=="N"))

coefs_for_plot_byposition <- bind_rows(
  tidy(ceda_fit_first_a) %>% mutate(subset = "2+ candidates"),
  tidy(ceda_fit_firstsecond_b) %>% mutate(subset = "3+ candidates"),
  tidy(ceda_fit_firstsecondthird_c) %>% mutate(subset = "4+ candidates"),
  tidy(ceda_fit_firstsecondthirdfourth_d) %>% mutate(subset = "5+ candidates")
) %>%
  filter(grepl("ballot_order", term)) %>%
  mutate(term_num = str_extract(term, "\\d+"),
         term_word = dplyr::recode(term_num,
                                   "1" = "First",
                                   "2" = "Second",
                                   "3" = "Third",
                                   "4" = "Fourth"),
         term_word = factor(term_word, levels=c("Fourth","Third","Second","First")),
         subset = factor(subset, levels=c("5+ candidates","4+ candidates","3+ candidates","2+ candidates"), ordered=T))

(figA8 <- ggplot(coefs_for_plot_byposition) +
  geom_hline(yintercept=0, lty=2, col="red") +
  geom_point(aes(x=term_word, y=estimate), position=position_dodge(width=1), size=2) +
  geom_errorbar(aes(y=estimate, x=term_word, ymin=estimate+qnorm(0.025)*std.error, ymax=estimate+qnorm(0.975)*std.error),
                position=position_dodge(width=1), width=0) +
  geom_errorbar(aes(y=estimate, x=term_word, ymin=estimate+qnorm(0.05)*std.error, ymax=estimate+qnorm(0.95)*std.error),
                lwd=1, width=0, position=position_dodge(width=1)) +
  facet_wrap(~subset, scales="free_y") +
  scale_y_continuous("Effect of being in ballot position on vote share", labels=scales::percent_format(accuracy = 1)) +
  scale_x_discrete("Ballot position") +
  coord_flip() +
  theme_minimal() +
  theme(text = element_text(family="Arial Narrow"),
        axis.title.y = element_text(margin = margin(r = 10), angle = 90),
        plot.margin = margin(.5, .5, 0, .5, "cm"),
        strip.background = element_rect(fill="grey"))
)

ggsave(figA8, filename = "figures/figureA8_effect_allpositions.png", height=6, width=6)

# ============================================================================
# Print completion message ===================================================
# ============================================================================
cat("\n\nReplication complete. All figures have been saved to the 'figures' subfolder.\n")
